home *** CD-ROM | disk | FTP | other *** search
/ Libris Britannia 4 / science library(b).zip / science library(b) / PROGRAMM / PASCAL / 0514.ZIP / CRAYZ15.ARC / DECOMP.FOR < prev    next >
Text File  |  1986-09-22  |  2KB  |  58 lines
  1. C COLLECTED ALGORITHMS FROM CACM
  2. C
  3. C ALGORITHM 423
  4. C LINEAR EQUATION SOLVER
  5. C CLEVE B. MOLER
  6. C DEPARTMENT OF MATHEMATICS, UNIVERSITY OF MICHIGAN
  7. C ANN ARBOR, MICHIGAN 48104
  8. C RECD. 1 JULY 1970 AND 1 DEC. 1970
  9. C
  10.       SUBROUTINE DECOMP(N, NDIM, A, IP)
  11.       REAL A(NDIM,NDIM),t
  12.       INTEGER IP(NDIM)
  13. C
  14. C  MATRIX TRIANGULARIZATION BY GAUSSIAN ELIMINATION.
  15. C  INPUT..
  16. C    N = ORDER OF MATRIX.
  17. C    NDIM = DECLARED DIMENSION OF ARRAY A .
  18. C    A = MATRIX TO BE TRIANGULARIZED.
  19. C  OUTPUT..
  20. C    A(I,J), I.LE.J = UPPER TRIANGULAR FACTOR, U .
  21. C    A(I,J), I.GT.J = MULTIPLIERS = LOWER TRIANGULAR
  22. C                     FACTOR, I-L .
  23. C    IP(K), K.LT.N = INDEX OF K-TH PIVOT ROW.
  24. C    IP(N) = (-1)**(NUMBER OF INTERCHANGES) OR 0 .
  25. C  USE 'SOLVE' TO OBTAIN SOLUTION OF LINEAR SYSTEM.
  26. C  DETERM(A) = IP(N)*A(1,1)*A(2,2)*...*A(N,N).
  27. C  IF IP(N)=0, A IS SINGULAR, SOLVE WILL DIVIDE BY ZERO.
  28. C  INTERCHANGES FINISHED IN U , ONLY PARTLY IN L .
  29. C
  30.       IP(N) = 1
  31.       DO 6 K = 1,N
  32.         IF(K.EQ.N) GO TO 5
  33.         KP1 = K+1
  34.         M = K
  35.         DO 1 I = KP1,N
  36.           IF(ABS(A(I,K)).GT.ABS(A(M,K))) M = I
  37.     1   CONTINUE
  38.         IP(K) = M
  39.         IF(M.NE.K) IP(N) = -IP(N)
  40.         T = A(M,K)
  41.         A(M,K) = A(K,K)
  42.         A(K,K) = T
  43.         IF(T.EQ.0.) GO TO 5
  44.         DO 2 I = KP1,N
  45.     2     A(I,K) = -A(I,K)/T
  46.         DO 4 J = KP1,N
  47.           T = A(M,J)
  48.           A(M,J) = A(K,J)
  49.           A(K,J) = T
  50.           IF(T.EQ.0.) GO TO 4
  51.           DO 3 I = KP1,N
  52.     3       A(I,J) = A(I,J) + A(I,K)*T
  53.     4   CONTINUE
  54.     5   IF(A(K,K).EQ.0.) IP(N) = 0
  55.     6 CONTINUE
  56.       RETURN
  57.       END
  58.